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Abstract 

We calculate numerically the probability e'x.p[jF{E / Esph, N/Ngph)] of the 
false vacuum decay in the massive four-dimensional —Xcp'^ model from multipar- 
ticle initial states with fixed number of particles A'^ and energy E greater than 
the height of the barrier Esph- We find that at E <^ SE'sp/i and N ^ OANgph 
the decay is classically forbidden and thus is exponentially suppressed. We 
argue that the classically forbidden region extends at small at least up to 
E ~ lOEsph and, most likely, to all energies. Our data suggest that the false 
vacuum decay induced by two-particle collisions is exponentially suppressed at 
all energies. 
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1 Introduction 



In the last few years a substantial progress has been made in understanding of non- 
vacuum tunneling in field theory. The interest to this subject is inspired by attempts 
to calculate, in the Standard Model, the instanton-mediated baryon number violation 
in particle collisions at high energy []I|. Although a variety of methods has been 
developed for calculating the tunneling probabihties at i? -C Egph, where Esph is the 
height of the barrier (the sphaleron energy in the Electroweak Theory) , the behaviour 
of probabilities in the most interesting region E > Egp^ is still unknown (for a review 
see ref.0). In particular, the question of whether the exponential suppression of 
baryon number violating processes disappears at some sufficiently high energy of 
colliding particles, still remains unanswered. Since analytical approaches to tunneling 
induced by particle collisions seem to be exhausted, in the present paper we address 
this problem by means of numerical methods. 

A suitable starting point for numerical study of non- vacuum tunneling is provided 
by the semiclassical approach developed in refs. 0, ^, |^ . In this approach the key role 
is played by the probability of tunneling from a mixed state with fixed energy E and 
number of particles A^. This probability is defined as follows, 

a{E,N)=J2\{f\SpEpN\^)\^ 

where S is the S-matrix, Pe,n are projectors onto subspaces of fixed energy E and 
fixed number of particles A^, respectively, while the states \i) and |/) are perturbative 
excitations above two vacua lying on different sides of the barrier. It was argued in 
refs.[^ ^ that at any fixed A^, the probability (y{E,N) sets an upper bound for the 
two-particle cross section, while in the limit of small A^ this probability reproduces 
the two-particle cross section with exponential accuracy. 

The advantage of considering multiparticle probability (t{E, N) instead of two- 
particle one is that in the weak coupling regime ^ and E,N ^ l/g"^, the former 
can be calculated semiclassically and has the form 

a(E,Ar)~exp{lF(e,z/)}, (1) 

where e = E/Egph, v = N/N^ph and Ngph ~ l/g'^ is the number of particles produced 
in the sphaleron decay. The function F{e, v) is negative at low energies, which corre- 
sponds to exponential suppression of probability in this domain. At e = 1 and v = 1 
(at the sphaleron), the function F(e, z/) vanishes, i.e. the exponential suppression 
disappears. The absence of exponential suppression of the two-particle cross section 
would show up as zero of the function F(e, v) at some fixed e > 1 and z/ ^ 0. 

The function F(e, u) is determined by a solution to the specific classical boundary 
value problem for the complexified field equations p. Namely, one has to solve the 
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usual field equation 

where time and the field are treated as complex variables. All the information about 
particular problem is encoded into the boundary conditions, which are formulated on 
the contour ABCD in the complex time plane (see Fig.l). The asymptotic regions A 
and D correspond to the initial and final states, respectively. 

The boundary conditions are as follows. The final boundary condition states that 
the field is real on the line CD, where it represents the classical evolution of the system 
after barrier penetration. This requirement can be imposed at the point C where it 
implies 

Im = 0, 

I.nf^O. ,3) 

The part CD of the contour is not essential for formulation of these boundary con- 
ditions. In fact, this part does not play any role in the calculation of the function 
F(e, u) and can be dropped, unless the details of the final state are of interest. Due 
to this fact the final boundary conditions apply also to — A0^ model where the 
stable vacuum does not exist and field develops a singularity in the final state, i.e. 
somewhere on the line CD (point P in Fig.l). 

At the other end of the contour, in the asymptotic region A, the field is required to 
be linear. The initial boundary conditions fix the ratio of amplitudes in the negative- 
and positive-frequency parts to be a constant independent of momentum. If one 
writes general asymptotics of the field in the form 



■I ^(2^0^ 



where = Ret and 6* is a real positive parameter, then the boundary condition is 

^7k = (/k)*. (5) 

This equation fixes one particular linear combination of the field and its time deriva- 
tive. In total, there are two real conditions at each end of the contour for the second 
order complex differential equation, so in general the boundary value problem is com- 
pletely specified. 

The input variables E and enter the above boundary conditions through two 
parameters: 9, which enters the boundary conditions explicitly, and T, which is 
the amount of Euclidean evolution (see Fig.l) and specifies the place where initial 
boundary conditions are imposed. The relation between T, 9 and E, N is given by 
the equations 

E = -^MiS{T,9)l 
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N=-2—Re[tSiT,e)], (6) 

where S{T, 9) is the action of the solution for given T and 6', evaluated along the con- 
tour. Note that according to boundary conditions the field is real on the line CD, 
so that this part of the contour does not contribute into Re[i5'(T, 6*)]. Alternatively, 
the energy and the number of particles can be read off from the initial asymptotics 
of the field, 

N = j dkf^f:, 

E = J dkujMt ■ (7) 

One can check that eqs.(|^) and coincide, provided that eqs.(||) and (|) are satisfied. 

Given the solution to the boundary value problem one can calculate the function 
-F(e, v) according to the formula 

^F{e,v) = 2ET + Ne + 2Re[iS{T,e)\, (8) 

where T and 9 depend on E and N through eqs.(P). 

Several remarks are in order. First, as follows from eqs.(|^) and (^, at 7^ 
the solution is necessarily complex along the line AB. Thus, contributions to the 
function F(e, v) come from both Minkowskian (AB) and Euclidean (BC) parts of 
the contour. Second, as the field equations are analytic in time, the solution is also 
analytic everywhere except possible singularities. Thus, the contour ABC can be 
deformed, provided that the asymptotic region A is untouched and singularities are 
not crossed. In fact, one must expect singularities in between the contour ABC and 
the negative part of the real time axis. Otherwise the reality of the field at real 
t and the initial boundary conditions are incompatible. In the model we consider 
below these singularities lie on the real time shown in Fig.l. In numerical 

calculations it is convenient to deform the contour in such a way that it passes far 
from singularities (dotted line in Fig.l). 

In one particular case, namely at = 0, the above boundary value problem 
simplifies considerably. In this case the initial boundary conditions, eq.(|), reduce to 
the reality condition imposed at Imt = T. The solution to the resulting boundary 
value problem is given by periodic instanton of ref.|P]. Periodic instanton is a real 
periodic solution to the Euclidean field equations with period 2T and two turning 
points at t = and t = iT. Being analytically continued to the Minkowskian domain 
through the turning points, periodic instanton stays real at lines Im t = and Im t = T 
and thus satisfies the boundary conditions with ^ = 0. 

The periodic instanton can be found analytically in two extreme cases: at e ^ 1 it 
can be approximated by the periodic chain of instantons and antiinstantons, while at 
1 — e <^ 1 it is approximately given by the oscillations in the sphaleron negative mode. 
At intermediate energies the periodic instanton can be obtained numerically |^ . In 
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the e-z/ plane, periodic instantons form a line connecting the points e = z/ = (zero 
energy instanton) and e = u = 1 (sphaleron). The function F(e, z/) monotonically 
grows from F(0, 0) = —Sinst to 1) = along this line. 

The rest of this paper is devoted to numerical solution of the boundary value 
problem specified by eqs.(^, (|^) and in the general case 6^ 7^ 0, and investigation 
of the behaviour of the function F(e, z/) at e ^ 1. For simplicity reasons we concentrate 
on particle-induced false vacuum decay in the theory with one scalar field. 



2 The Model 

The choice of particular model turns out to be strongly constrained for technical 
reasons. After several attempts we ended up with four- dimensional —Xcj)'^ theory 
with the mass term. The action of this model reads 

S = j dS{\d,<pd^<P - \m^^ + (9) 

where A is a positive constant. At non-zero m the region and the instability 
region > mj^JX are separated by a finite energy barrier, and one can ask whether 
the presence of colliding particles enhances the probability of the barrier penetration 
(i.e., the decay of the metastable state = 0). 

It is a lucky coincidence that this model is of special interest for another reason, 
as it has many features reminiscent of bosonic sector of the Electroweak Theory. At 
m = (zero Higgs vacuum expectation value) both models are conformally invariant 
and possess instanton solutions. In the case of — A0^ model the instanton ||^ has the 
action ^ 

Sinst = (10) 

and describes the decay of the metastable state = 0. At m 7^ (non-zero Higgs 
vacuum expectation value) the conformal symmetry is softly broken and instanton 
solutions disappear, while the low energy transitions are described by constrained 
instantons [|10[- The action ([T0|) gets small energy-dependent correction and still de- 
termines the transition probability. Like the Electroweak Theory, the — A0^ model 
possesses the sphaleron solution. The energy of the sphaleron can be found numeri- 
cally, 

TTl 

Esph = 18.9y. (11) 

Before discussing the discrete formulation of the boundary value problem specified 
above, it is convenient to rewrite the action (^) in the dimensionless variables. Since 
the problem is 0(3)-symmetric, we restrict ourselves to s-wave scattering. The change 
of variables 

1 



\x 
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X = m ^y, 

brings the action for spherically symmetric configurations to the form 

no 

1 o 1 



JdtJ dr[-{d,ijr - ^{dri^r - ^i^' + (12) 



where r = \y\ is the dimensionless radial distance. Throughout the rest of this paper 
all dimensionfull quantities are measured in the units of mass. 

The starting point for our calculations is the discretized version of the action 
(|T2|). The system is put on a grid of the size L with Ux + I sites rj = jL/rix, 
j = . . .Ux- Similarly, the time contour is represented by the set of complex points 
tj with i = . . .nt- The field ip{t, r) transforms into the set of complex variables ipij- 
To define integrations, we introduce two sets of intervals for each coordinate, 

drj = Tj+i -rj, j = 0...nx- 1, 

drj = {drj^i + drj)/2, j = 1 . . .n^ - I, 
dr o,n, = dro,nj2, 

and similarly for dti and dti (tilted and non-tilted intervals are used to integrate fields 
and derivatives, respectively). With these definitions, the discretized action reads 

= y - i^ij) - 2 (^^.i+i - ^^i) ^ - Vijdtidrj\ , (13) 

ij ^ 3 

where ^ ^ 

^ii = lfl-^'Ptj at jVO, 

^.0 = 0. 

Equations of motion can be easily derived from this action by variation with respect 

to 

We also need the discrete version of the boundary conditions. Let us start with 
the final boundary conditions. The first equation (|^) transforms into 

Im-ipntj = . 

The discrete analog of the second equation can be derived from the requirement 
that the solution, being continued to the real time axis according to the discretized 
equations of motion, stays real. This requirement leads to the equations 

T n 

9lpnt,j 
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The discrete reformulation of the initial boundary conditions requires somewhat 
more work. It can be obtained by means of the following trick. Let us make one step 
back and treat time as continuous variable. Then we would obtain the n^-dimensional 



quantum mechanical system with the action trivially derived from eq . (pTSf) . Since it 
is assumed that the field reaches linear regime in the initial asymptotic region, we 
can restrict ourselves to the quadratic part of the resulting action and diagonalize it 
numerically. In this way we obtain the discrete analog of the momentum represen- 
tation. Given the corresponding eigenvalues and eigenvectors, it is straightforward 
to perform the decomposition of the field into positive- and negative-frequency parts 
and impose boundary conditions (|^). Eqs.(|^) are then translated into conditions 
imposed on some linear combinations of tpoj and tpij. These are the desired initial 
boundary conditions in the discrete formulation. To save space we do not present the 
corresponding cumbersome expressions. 

It is instructive to note that the above n^j-dimensional quantum mechanical system 
itself can be viewed as a model for studying tunneling transitions from the excited 
states in multidimensional systems. The action (|T3p can be treated as discretized 
(in time) version of the action for this quantum mechanical system. In practice this 
means that particular number of lattice sites in space direction is not crucial for our 
conclusions, as long as linear regime can be reached at the initial part of the contour. 



3 Numerical method 

To solve the boundary value problem for given values of T and 6 we use a multidi- 
mensional version of Newton's method. This is a relaxation procedure which takes 
as input an approximate solution to the non-linear field equations and improves it at 
each iteration by solving the linearized equations in the background of the current 
approximation. Iterations are repeated until non-linear equations are satisfied to the 
desired accuracy. The advantage of this algorithm is that its convergency requires 
neither positive-definiteness nor even reality of the action. Moreover, the convergency 
is quadratic, provided that initial approximation is choosen sufficiently close to the 
solution. In practice, the accuracy of 10~^° is reached in 3-6 iterations. 

A drawback of the Newton's method is that its basin of convergency can be very 
narrow. Thus, it is important to have a good initial approximation at least for some 
values of T and 6. Then one can move gradually in the T - 9 plane using output of 
each run as input for the next one. In our case the starting configuration is provided 
by the periodic instanton which, aX E Egph, can be approximated by sphaleron 
plus harmonic oscillation in the sphaleron negative mode (both the sphaleron and its 
negative mode have to be found numerically). The period of the periodic instanton 
at E —>■ Esph is determined by the sphaleron negative eigenvalue and equals 

Tcrit = 0.78 . 

The amplitude of the oscillation goes to zero when T Tcrit- So, we take T ^ Tcrit, 
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^ = 0, and adjust the amplitude of the oscillation in order to get inside the basin of 
convergency of the periodic instanton with the period T. In this way we obtain the 
very first solution to our boundary value problem. Changing T by small steps we then 
reproduce the whole family of periodic instantons with different periods. Calculations 
show that in the — model the period of the periodic instanton varies from Tcru to 
zero when the energy decreases from Egph to zero. This unusual behaviour of period 
with energy is much in common with that in the 0(3) sigma model [Q, where it is 
also related to softly broken conformal invariance. On the basis of this similarity one 
should expect the analogous behaviour in the Electroweak Theory. 

The Newton's method reduces the non-linear boundary value problem to sequen- 
tial solution of a few linear boundary value problems. At each iteration, to find the 
correction 6ip to the current approximation amounts to one matrix inversion, 

6ip = D-^R , 

where D is a matrix of second derivatives of the action and R is a. vector of first 
derivatives, both evaluated at the background. The matrix D has the dimension 
{2ntnx) X [2ntnx), the factor 2 being due to the complexity of the field. In general, 
the computer time necessary for inversion of such a matrix scales like (n(?T,^)^. One 
can use, however, the fact that this matrix originates from the local second order 
differential equation and thus is sparse. The sparseness enables one to invert this 
matrix in oc ntn^ operations. Note that for systems containing nf fields instead of 
one, this number would be nt{nfnx)^. In practice, to find one solution to the boundary 
value problem at the grid of dimension rit x = 200 x 40 takes about 5 minutes at 
SPARCstation 20, most of this time being spent for matrix inversion. The necessary 
amount of memory scales like ritnl and for the above grid is of order 10 Mb. 

Let us now discuss constraints on other grid parameters. The most restrictive 
requirement is that the field must reach linear regime in the asymptotic region A. 
For that the number of independent modes (which equals ra^) must be large enough, 
and the length of Minkowskian part of the contour, Tm, must be sufficient to let the 
energy spread over these modes. Since in the Minkowskian region the most of the 
field propagates along the light cone, the space size L should be larger than Tm to 
avoid boundary effects. On the other hand, the lattice spacing in r direction, Ar, 
should be small enough not only to be close to the continuum limit, but also because 
otherwise the free spectrum would be cut at too low frequency, ujmax ~ vr/Ar. This 
would impose constraints on the available region in the E-N space, E/N <^ i^max, 
and would not allow for the initial states consisting of small number of high energy 
particles. Thus, all these arguments push towards large n^, which is however bound 
at relatively low value n^; ~ 50 — 100 by the abilities of available computers. 

The number of time slices rit is not so constrained since the required computer 
resources scale linearly with rij. In our calculations we take rif few times larger than 
in order to saturate the continuum limit in time. This number is different for different 
T and 6. 
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The particular choice of parameters is model-dependent and is made by trial 
and error. In the model we consider here rather small value Tm = 3 is sufficient 
to reach linear regime starting from the vicinity of the sphaleron configuration. Fast 
linearization is the advantage of high-dimensional models [jll| and is due to the volume 



factor (note in the interaction term in eg. ([T^) ) . This effect is absent in two- 
dimensional models. Moreover, in — model, the value Tm = 3 is sufficient for 
linearization at least up to energies E ~ SEgph- This feature differs the above model 
from other four-dimensional models we have tried. 

Fast linearization allows to take relatively small space size, L = 3, which leads to 
the sufficiently wide spectrum already at n^, = 40. The free spectrum for these values 
of parameters reproduces the continuum spectrum with reasonable accuracy up to 
u ~ 15. 



4 Results and discussion 

The results presented in this this paper were obtained at L = 3, Tm = 3, = 40 
and rit varying from 200 to 300. For these values of grid parameters we have found 
about 1500 solutions to the boundary value problem with different E and A^. The 
region of e-u plane covered by the solutions is shown in Fig.2. The top left corner of 
the solution region (point S in Fig.2) is the sphaleron. It corresponds to e = u = 1. 
At this point the function F(e, u) vanishes and exponential suppression of probability 
disappears. 

Towards small energies the solution region is bound by the line of periodic in- 
stantons (the line SP in Fig.2). In the continuum limit this line would end up at 
the point e = z/ = which corresponds to the zero energy instanton. However, the 
spatial size of periodic instantons rapidly decreases along this line. At the above val- 
ues of grid parameters, the instanton size becomes comparable to the lattice spacing 
at e ~ 0.4 (the point P in Fig.2). At this point the Newton's algorithm stops to 
converge. Along the line of periodic instantons, the function F{e, v) monotonically 
decreases and reaches the value F = —O.GSinst at the point P. 

The boundary of convergency continues at approximately constant u towards 
higher energy. It is represented by the line PQ in Fig.2. The function F is nega- 
tive along this line. 

The third boundary of the solution region is formed by the line F = (solid line in 
Fig.2). It starts at the sphaleron and goes towards higher energy and smaller number 
of particles. At energy E ^ SEgph it comes close to the boundary of convergency, 
which makes finite the available part of e-u plane. Our data cover this available 
region. 

In this paper we concentrate on the dependence of the function F on energy and 
number of particles. It is most conveniently represented by the lines of constant F in 
the e-u plane. These lines can be obtained by interpolation and are shown in Fig.3. 
They start at the line of periodic instantons with infinite negative slope. The latter 
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can be seen analytically. Indeed, from eqs. (|^) and (§) one immediately derives that 



dN 
'dE F 



2T 



(14) 



The slope dN/dE\p becomes infinite at the periodic instantons where ^ = 0. 

The function F(e, v) determines the maximum probability of induced false vac- 
uum decay among A^-particle initial states. The line F = separates the classically 
forbidden and classically allowed regions. We see from Fig. 3 that for each energy 
E > Esph there exists a minimum number of particles Ncrit{E) for which the decay 
may occur classically without exponential suppression. At > Ncru one thus may 
expect the existence of classical configurations which pass above the barrier. On the 
other hand, at iV < Ncru all classical solutions bounce off the barrier. By continuity, 
when approaches Ncru from above, the classical solutions which pass above the 
barrier should spend more and more time on its top oscillating above the sphaleron, 
so that in the limit A^ = Ncru this time goes to infinity. 

At A^ < Ncrit, the decay is a tunneling event which is described by corresponding 
solution to our boundary value problem. As discussed above, the real time part of 
this tunneling solution represents the classical evolution of the system after barrier 
penetration. In our model this evolution leads to the singularity (point P in Fig.l). In 
between two singularities on the real time axis, the tunneling solution represents the 
field which comes from infinity, bounces off the barrier and goes back to infinity. When 
A^ approaches Ncru from below, the distance between two singularities grows (see 
Fig. 4) and the solution spends more and more time around the sphaleron. Although 
due to instabilities we were not able to trace the positions of singularities when the 
distance between them becomes of order one, we expect that at the line F = this 
distance goes to infinity and the real time part of the tunneling solution reproduces 
the infinitely oscillating classical solution described above. Thus, the configurations 
which correspond to the line F = should be classical solutions infinitely oscillating 
above the sphaleron. 

Due to the existence of minimum number of particles at fixed energy E, one can 
try to obtain the line F = by minimizing A^ over the set of classical solutions 
which pass above the barrier. This strategy was used in ref . [0 to set bounds on the 
forbidden region. We have tried to use this strategy to reproduce the line F = found 
by the tunneling approach. We have obtained bounds which are noticeably higher 
than the real position of this line. The reason could be that the classical solutions 
which correspond to the line F = are infinitely oscillating around the sphaleron and 
thus are singular points of the space of solutions which pass above the barrier. 

As one can see from Fig. 3, the size of allowed region increases towards higher en- 
ergy, while the minimum number of particles for which the decay is not exponentially 
suppressed, decreases. We have traced this behaviour up to F ~ 3i?sp/i where the 
absence of convergency prevents us from going further with the same grid parameters. 
Since the probability (j{E, N) is an upper bound for the two-particle cross section at 
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energy E, we conclude that the latter cross section is exponentially suppressed at 
E < 3Esph. Moreover, extrapolating the behaviour of NcritiE) and assuming that 
the derivative dNcrit{E) / dE continues to decrease in absolute value, one can see that 
the line F = does not cross the = axis at least up to ~ lOEgph- Thus, 
we conclude that exponential suppression of the two-particle cross section persists at 
least up to ii^ ~ lOEgph- 

Careful look at the data shows that the behaviour of the function F{e, v) may 
qualitatively change at substantially lower energy. To see this consider the slope of 
the lines of constant F as a function of E. The derivative dN/dE\p can be expressed 
through the known parameters according to eq.(|l^). Its behaviour with energy is 
shown in Fig. 5. Extrapolating, one can see that it would become zero at E^ ~ t^.SEgpfi 
and A^* ~ 0.35Nsph- If this indeed happens, then starting at this energy the line F = 
in coordinates E and would turn up {dNcrit{E) / dE would become positive), while 
the function F{E, N^,) aX E > E^ would decrease with energy. Clearly, this behaviour 
is unphysical. Indeed, it is always possible to put the excess energy E — E^ m. one 
quantum particle and thus effectively decrease the energy involved in the semiclassical 
tunneling, leaving the number of particles practically unchanged. The corresponding 
initial state would not be semiclassical, and its probability of tunneling would be 
F{E^, N^:), i.e. higher than follows from the semiclassical picture. If this situation 
realizes, it is natural to expect that starting at E = E^ the correct function F{E, N^) 
stays constant and thus the line F = never crosses the N = axis. The latter 
would mean that the two-particle cross section is always exponentially suppressed. 

As is clear from Fig.3, in order to go to higher energies and verify the above 
conclusions based on extrapolation, one has to extend the convergency region towards 
smaller number of particles. The breakdown of convergency is most likely related to 
the lattice spacing in r direction. We have checked that doubling at fixed L 
extends the convergency region for periodic instanton solutions. We expect that the 
same effect persists at 6' 7^ 0. Thus, increasing Ux could be the way to extend the 
convergency region and improve our results. This, however, will require noticeably 
more computer resources. 

In conclusion, we would like to stress again the analogy, as long as instanton-like 
transitions are concerned, between the — model and the bosonic sector of the 
Electroweak Theory. On the basis of this analogy, we expect that the Electroweak 
Theory is also a suitable model for studying instanton-like transitions numerically. 
Since, because of larger number of fields, calculations there are at least 4^ times slower, 
they would require the use of computers much faster than a typical workstation. 
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Figure captions 



Fig.l. The contour ABCD in the complex time plane where the boundary value 
problem is formulated. Crossed circles represent singularities of the field. Dotted line 
schematically shows the deformed contour used in numerical calculations. 

Fig.2. The e-v map of the obtained solutions. The point S corresponds to the 
sphaleron. Crossed circles lying on the line SP are periodic instantons. The line PQ 
represents the boundary of convergency region, while the solid line corresponds to 
F(e,z/) = 0. 

Fig.3. The hues of constant F(e, v) in the e-v plane. Numbers show the value of F 
in the units of Sinst- 

Fig. 4. The distance between singularities At as a function of F when F approaches 
zero. Dashed line corresponds to the region where the procedure of calculating At is 
not rehable. 

Fig. 5. The dependence of —dN/dE\p = 2T/9 on e. Dashed line represents the 
extrapolation to higher energies, = E^/Egph corresponds to the point where the 
behaviour of F(e, v) may qualitatively change. 
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Fig.5 
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